home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libblas / src_original / zhbmv.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  10.0 KB  |  313 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE ZHBMV ( UPLO, N, K, ALPHA, A, LDA, X, INCX,
  5.      $                   BETA, Y, INCY )
  6. *     .. Scalar Arguments ..
  7.       COMPLEX*16         ALPHA, BETA
  8.       INTEGER            INCX, INCY, K, LDA, N
  9.       CHARACTER*1        UPLO
  10. *     .. Array Arguments ..
  11.       COMPLEX*16         A( LDA, * ), X( * ), Y( * )
  12. *     ..
  13. *
  14. *  Purpose
  15. *  =======
  16. *
  17. *  ZHBMV  performs the matrix-vector  operation
  18. *
  19. *     y := alpha*A*x + beta*y,
  20. *
  21. *  where alpha and beta are scalars, x and y are n element vectors and
  22. *  A is an n by n hermitian band matrix, with k super-diagonals.
  23. *
  24. *  Parameters
  25. *  ==========
  26. *
  27. *  UPLO   - CHARACTER*1.
  28. *           On entry, UPLO specifies whether the upper or lower
  29. *           triangular part of the band matrix A is being supplied as
  30. *           follows:
  31. *
  32. *              UPLO = 'U' or 'u'   The upper triangular part of A is
  33. *                                  being supplied.
  34. *
  35. *              UPLO = 'L' or 'l'   The lower triangular part of A is
  36. *                                  being supplied.
  37. *
  38. *           Unchanged on exit.
  39. *
  40. *  N      - INTEGER.
  41. *           On entry, N specifies the order of the matrix A.
  42. *           N must be at least zero.
  43. *           Unchanged on exit.
  44. *
  45. *  K      - INTEGER.
  46. *           On entry, K specifies the number of super-diagonals of the
  47. *           matrix A. K must satisfy  0 .le. K.
  48. *           Unchanged on exit.
  49. *
  50. *  ALPHA  - COMPLEX*16      .
  51. *           On entry, ALPHA specifies the scalar alpha.
  52. *           Unchanged on exit.
  53. *
  54. *  A      - COMPLEX*16       array of DIMENSION ( LDA, n ).
  55. *           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
  56. *           by n part of the array A must contain the upper triangular
  57. *           band part of the hermitian matrix, supplied column by
  58. *           column, with the leading diagonal of the matrix in row
  59. *           ( k + 1 ) of the array, the first super-diagonal starting at
  60. *           position 2 in row k, and so on. The top left k by k triangle
  61. *           of the array A is not referenced.
  62. *           The following program segment will transfer the upper
  63. *           triangular part of a hermitian band matrix from conventional
  64. *           full matrix storage to band storage:
  65. *
  66. *                 DO 20, J = 1, N
  67. *                    M = K + 1 - J
  68. *                    DO 10, I = MAX( 1, J - K ), J
  69. *                       A( M + I, J ) = matrix( I, J )
  70. *              10    CONTINUE
  71. *              20 CONTINUE
  72. *
  73. *           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
  74. *           by n part of the array A must contain the lower triangular
  75. *           band part of the hermitian matrix, supplied column by
  76. *           column, with the leading diagonal of the matrix in row 1 of
  77. *           the array, the first sub-diagonal starting at position 1 in
  78. *           row 2, and so on. The bottom right k by k triangle of the
  79. *           array A is not referenced.
  80. *           The following program segment will transfer the lower
  81. *           triangular part of a hermitian band matrix from conventional
  82. *           full matrix storage to band storage:
  83. *
  84. *                 DO 20, J = 1, N
  85. *                    M = 1 - J
  86. *                    DO 10, I = J, MIN( N, J + K )
  87. *                       A( M + I, J ) = matrix( I, J )
  88. *              10    CONTINUE
  89. *              20 CONTINUE
  90. *
  91. *           Note that the imaginary parts of the diagonal elements need
  92. *           not be set and are assumed to be zero.
  93. *           Unchanged on exit.
  94. *
  95. *  LDA    - INTEGER.
  96. *           On entry, LDA specifies the first dimension of A as declared
  97. *           in the calling (sub) program. LDA must be at least
  98. *           ( k + 1 ).
  99. *           Unchanged on exit.
  100. *
  101. *  X      - COMPLEX*16       array of DIMENSION at least
  102. *           ( 1 + ( n - 1 )*abs( INCX ) ).
  103. *           Before entry, the incremented array X must contain the
  104. *           vector x.
  105. *           Unchanged on exit.
  106. *
  107. *  INCX   - INTEGER.
  108. *           On entry, INCX specifies the increment for the elements of
  109. *           X. INCX must not be zero.
  110. *           Unchanged on exit.
  111. *
  112. *  BETA   - COMPLEX*16      .
  113. *           On entry, BETA specifies the scalar beta.
  114. *           Unchanged on exit.
  115. *
  116. *  Y      - COMPLEX*16       array of DIMENSION at least
  117. *           ( 1 + ( n - 1 )*abs( INCY ) ).
  118. *           Before entry, the incremented array Y must contain the
  119. *           vector y. On exit, Y is overwritten by the updated vector y.
  120. *
  121. *  INCY   - INTEGER.
  122. *           On entry, INCY specifies the increment for the elements of
  123. *           Y. INCY must not be zero.
  124. *           Unchanged on exit.
  125. *
  126. *
  127. *  Level 2 Blas routine.
  128. *
  129. *  -- Written on 22-October-1986.
  130. *     Jack Dongarra, Argonne National Lab.
  131. *     Jeremy Du Croz, Nag Central Office.
  132. *     Sven Hammarling, Nag Central Office.
  133. *     Richard Hanson, Sandia National Labs.
  134. *
  135. *
  136. *     .. Parameters ..
  137.       COMPLEX*16         ONE
  138.       PARAMETER        ( ONE  = ( 1.0D+0, 0.0D+0 ) )
  139.       COMPLEX*16         ZERO
  140.       PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
  141. *     .. Local Scalars ..
  142.       COMPLEX*16         TEMP1, TEMP2
  143.       INTEGER            I, INFO, IX, IY, J, JX, JY, KPLUS1, KX, KY, L
  144. *     .. External Functions ..
  145.       LOGICAL            LSAME
  146.       EXTERNAL           LSAME
  147. *     .. External Subroutines ..
  148.       EXTERNAL           XERBLA
  149. *     .. Intrinsic Functions ..
  150.       INTRINSIC          DCONJG, MAX, MIN, DBLE
  151. *     ..
  152. *     .. Executable Statements ..
  153. *
  154. *     Test the input parameters.
  155. *
  156.       INFO = 0
  157.       IF     ( .NOT.LSAME( UPLO, 'U' ).AND.
  158.      $         .NOT.LSAME( UPLO, 'L' )      )THEN
  159.          INFO = 1
  160.       ELSE IF( N.LT.0 )THEN
  161.          INFO = 2
  162.       ELSE IF( K.LT.0 )THEN
  163.          INFO = 3
  164.       ELSE IF( LDA.LT.( K + 1 ) )THEN
  165.          INFO = 6
  166.       ELSE IF( INCX.EQ.0 )THEN
  167.          INFO = 8
  168.       ELSE IF( INCY.EQ.0 )THEN
  169.          INFO = 11
  170.       END IF
  171.       IF( INFO.NE.0 )THEN
  172.          CALL XERBLA( 'ZHBMV ', INFO )
  173.          RETURN
  174.       END IF
  175. *
  176. *     Quick return if possible.
  177. *
  178.       IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
  179.      $   RETURN
  180. *
  181. *     Set up the start points in  X  and  Y.
  182. *
  183.       IF( INCX.GT.0 )THEN
  184.          KX = 1
  185.       ELSE
  186.          KX = 1 - ( N - 1 )*INCX
  187.       END IF
  188.       IF( INCY.GT.0 )THEN
  189.          KY = 1
  190.       ELSE
  191.          KY = 1 - ( N - 1 )*INCY
  192.       END IF
  193. *
  194. *     Start the operations. In this version the elements of the array A
  195. *     are accessed sequentially with one pass through A.
  196. *
  197. *     First form  y := beta*y.
  198. *
  199.       IF( BETA.NE.ONE )THEN
  200.          IF( INCY.EQ.1 )THEN
  201.             IF( BETA.EQ.ZERO )THEN
  202.                DO 10, I = 1, N
  203.                   Y( I ) = ZERO
  204.    10          CONTINUE
  205.             ELSE
  206.                DO 20, I = 1, N
  207.                   Y( I ) = BETA*Y( I )
  208.    20          CONTINUE
  209.             END IF
  210.          ELSE
  211.             IY = KY
  212.             IF( BETA.EQ.ZERO )THEN
  213.                DO 30, I = 1, N
  214.                   Y( IY ) = ZERO
  215.                   IY      = IY   + INCY
  216.    30          CONTINUE
  217.             ELSE
  218.                DO 40, I = 1, N
  219.                   Y( IY ) = BETA*Y( IY )
  220.                   IY      = IY           + INCY
  221.    40          CONTINUE
  222.             END IF
  223.          END IF
  224.       END IF
  225.       IF( ALPHA.EQ.ZERO )
  226.      $   RETURN
  227.       IF( LSAME( UPLO, 'U' ) )THEN
  228. *
  229. *        Form  y  when upper triangle of A is stored.
  230. *
  231.          KPLUS1 = K + 1
  232.          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  233.             DO 60, J = 1, N
  234.                TEMP1 = ALPHA*X( J )
  235.                TEMP2 = ZERO
  236.                L     = KPLUS1 - J
  237.                DO 50, I = MAX( 1, J - K ), J - 1
  238.                   Y( I ) = Y( I ) + TEMP1*A( L + I, J )
  239.                   TEMP2  = TEMP2  + DCONJG( A( L + I, J ) )*X( I )
  240.    50          CONTINUE
  241.                Y( J ) = Y( J ) + TEMP1*DBLE( A( KPLUS1, J ) )
  242.      $                         + ALPHA*TEMP2
  243.    60       CONTINUE
  244.          ELSE
  245.             JX = KX
  246.             JY = KY
  247.             DO 80, J = 1, N
  248.                TEMP1 = ALPHA*X( JX )
  249.                TEMP2 = ZERO
  250.                IX    = KX
  251.                IY    = KY
  252.                L     = KPLUS1 - J
  253.                DO 70, I = MAX( 1, J - K ), J - 1
  254.                   Y( IY ) = Y( IY ) + TEMP1*A( L + I, J )
  255.                   TEMP2   = TEMP2   + DCONJG( A( L + I, J ) )*X( IX )
  256.                   IX      = IX      + INCX
  257.                   IY      = IY      + INCY
  258.    70          CONTINUE
  259.                Y( JY ) = Y( JY ) + TEMP1*DBLE( A( KPLUS1, J ) )
  260.      $                           + ALPHA*TEMP2
  261.                JX      = JX      + INCX
  262.                JY      = JY      + INCY
  263.                IF( J.GT.K )THEN
  264.                   KX = KX + INCX
  265.                   KY = KY + INCY
  266.                END IF
  267.    80       CONTINUE
  268.          END IF
  269.       ELSE
  270. *
  271. *        Form  y  when lower triangle of A is stored.
  272. *
  273.          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  274.             DO 100, J = 1, N
  275.                TEMP1  = ALPHA*X( J )
  276.                TEMP2  = ZERO
  277.                Y( J ) = Y( J ) + TEMP1*DBLE( A( 1, J ) )
  278.                L      = 1      - J
  279.                DO 90, I = J + 1, MIN( N, J + K )
  280.                   Y( I ) = Y( I ) + TEMP1*A( L + I, J )
  281.                   TEMP2  = TEMP2  + DCONJG( A( L + I, J ) )*X( I )
  282.    90          CONTINUE
  283.                Y( J ) = Y( J ) + ALPHA*TEMP2
  284.   100       CONTINUE
  285.          ELSE
  286.             JX = KX
  287.             JY = KY
  288.             DO 120, J = 1, N
  289.                TEMP1   = ALPHA*X( JX )
  290.                TEMP2   = ZERO
  291.                Y( JY ) = Y( JY ) + TEMP1*DBLE( A( 1, J ) )
  292.                L       = 1       - J
  293.                IX      = JX
  294.                IY      = JY
  295.                DO 110, I = J + 1, MIN( N, J + K )
  296.                   IX      = IX      + INCX
  297.                   IY      = IY      + INCY
  298.                   Y( IY ) = Y( IY ) + TEMP1*A( L + I, J )
  299.                   TEMP2   = TEMP2   + DCONJG( A( L + I, J ) )*X( IX )
  300.   110          CONTINUE
  301.                Y( JY ) = Y( JY ) + ALPHA*TEMP2
  302.                JX      = JX      + INCX
  303.                JY      = JY      + INCY
  304.   120       CONTINUE
  305.          END IF
  306.       END IF
  307. *
  308.       RETURN
  309. *
  310. *     End of ZHBMV .
  311. *
  312.       END
  313.